function [Tinert2crtbp,Tcrtbp2inert,S_crtbp] = Inertial2CRTBP(JD_interest,S_eci,JD,R2,V2,body)

%% Constants
G           = 6.674e-11 * ((1/1000)^3);
mass_Mars   = 6.39e23;
mass_Deimos = 1.8e15;
mass_Phobos = 10.8e15;
GM_Mars     = G*mass_Mars;
GM_Deimos   = G*mass_Deimos;
GM_Phobos   = G*mass_Phobos;
GM_sun      = 1.327124400*10^(11);
%[km^3/s^2] Gravitational Parameters

u_sm = GM_Mars / (GM_Mars+GM_sun);
u_md = GM_Deimos / (GM_Deimos+GM_Mars);
u_mp = GM_Phobos / (GM_Phobos+GM_Mars);
%[] Mu constant
%    sm --> Sun-Mars system
%    em --> Earth-Moon system%%

DU_sm = 235.34e6;
DU_mp = (9236.53+9517.58)/2;
DU_md = (23455.5+23470.9)/2;
%[km] Distance from Sun to Mars

TU_sm = sqrt(DU_sm^3/GM_sun);
TU_mp = sqrt(DU_mp^3/GM_Mars);
TU_md = sqrt(DU_md^3/GM_Mars);
%[s] Time unit
%[] Parameter for altering the halo orbit's state
%%

if body == 'Phobos'
    
    GM_1= GM_Mars;          %Earth
    GM_2 = GM_Phobos;       %Phobos
    TU = TU_mp;
    LU = DU_mp;
    VU = LU/TU;
    Wcm = sqrt((GM_1+GM_2)/LU^3)*[0;0;1];
    
elseif body == 'Deimos'
    GM_1= GM_Mars;          %Mars
    GM_2 = GM_Deimos;       %Deimos
    TU = TU_md;             % Time unit
    LU = DU_md;             % Distance Unit
    VU = LU/TU;             % Velocity Unit
    Wcm = sqrt((GM_1+GM_2)/LU^3)*[0;0;1]; % Angular velocity of Deimos in MCI frame.
end
%% 1: Using julian date, determine moon's state at interested julian date.

JD_difference_vector = JD-JD_interest;
%[] Determine the difference time vector

ind = find(JD_difference_vector > 0 , 1 , 'first') - 1;
%[] Find the index of most close previous state

R2_o = R2(:,ind);
V2_o = V2(:,ind);
S2_o = [R2_o;V2_o];
%[km,km/s] Position and velocity of the moon WRT Earth.

UTC_JD_interest = CalendarDate(JD_interest);
UTC_JD_closest = CalendarDate(JD(ind));
%[] Get UTC time of JD time of interest, which is the JD_interest and
%JD(ind), which is the closest previous information.

% time_difference_seconds = etime(UTC_JD_closest,UTC_JD_interest);
time_difference_seconds = etime(UTC_JD_interest,UTC_JD_closest);
%[] Determine time difference in seconds.

if time_difference_seconds == 0
    R2_interest = R2_o;
    V2_interest = V2_o;
    S2_interest = [R2_interest;V2_interest];
    %[] Position and velocity of the moon at JD_interest.
else
    
    [~,S] = ode45(@(t,S)R2bpCartesianModel(t,S,GM_1),[0,time_difference_seconds],S2_o);
    S = S';
    %[] Integrate for difference amount to determine exact state of the moon
    %wrt Earth in ECI frame.
    
    R2_interest = S(1:3,end);
    V2_interest = S(4:6,end);
    S2_interest = [R2_interest;V2_interest];
    %[] Position and velocity of the moon at JD_interest.
end
%% 2: Determine states of CM

Scm1 = (GM_2*S2_interest)/(GM_2+GM_1);
Rcm1 = Scm1(1:3);
% Vcme = Scme(4:6);
%[] Determine the state of the center of mass at JD of interest
%There are no additional terms because the Earth's position wrt Earth is
%zero vector and the gravitational parameter of the satellite is assumed to
%be zero.

%% 3: Make all states WRT CM

% S1cm = -Scm1;
% Rearthcm = Searthcm(1:3);
% Vearthcm = Searthcm(4:6);
%[] State of the Earth wrt CM

% S2cm = S2_interest-Scm1;
% Rmooncm = Smooncm(1:3);
% Vmooncm = Smooncm(4:6);
%[] State of the moon wrt CM

Ssatcm = S_eci-Scm1;
Rsatcm = Ssatcm(1:3);
Vsatcm = Ssatcm(4:6);
%[] State of the state of interest [S_eci] wrt Center of mass.

%% 4: Determine Teci2crtbp and Tcrtbp2eci

Hm = cross(R2_interest,V2_interest);
%[1/s]Current specific angular momentum of the Moon WRT the CM in ECI coordinates.

Rhat = R2_interest / norm(R2_interest);
%[]Current CRTBP radial direction in ECI coordinates.

Nhat = Hm / norm(Hm);
%[]Current CRTBP normal direction in ECI coordinates.

That = cross(Nhat,Rhat);
%[]Current CRTBP tangential direction in ECI coordinates.

Tcrtbp2inert = [Rhat, That, Nhat];
%[]Matrix that transforms vectors from CRTBP to ECI coordinates.

Tinert2crtbp = transpose(Tcrtbp2inert);
%[] Matrix that transforms vectors from ECI to CRTBP coordinates.


%% 5: Determine S_crtbp

R_crtbp = Tinert2crtbp * Rsatcm;
%[] Determine positon of interest in CRTBP frame.

V_crtbp = Tinert2crtbp * Vsatcm - cross(Wcm,Tinert2crtbp*Rcm1) - cross(Wcm,R_crtbp);
%[] Determine positon of interest in CRTBP frame.

S_crtbp = [R_crtbp/LU;V_crtbp/VU];
%[] State of the S_eci eexpressed at CRTBP frame at JD of interest.

end